Pollutant inventory spatialization

1A4-Residential-Tertiary

This document provides the methodlogy and the main results regarding the spatialization of pollutation inventory.

1A4bi - Residential: Stationary combustion

source.1A4bi <- list(sources = list(points = NA, lines = NA, polygon = NA), total = list(spatialize = NA, inventory = NA))

source.1A4bi$total$spatialize <- readxl::read_xlsx(path = source.file, range = "D32:I32", sheet = source.sheet, col_names = vars)
source.1A4bi$total$inventory <- readxl::read_xlsx(path = source.file, range = "D37:I37", sheet = source.sheet, col_names = vars)

clc_18 <- readOGR("D:/R_projects/Spatialization/Data/clc/CLC18_RS.shp")

## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS, : Discarded datum European_Terrestrial_Reference_System_1989 in CRS
## definition: +proj=utm +zone=34 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
## +no_defs

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\R_projects\Spatialization\Data\clc\CLC18_RS.shp", layer: "CLC18_RS"
## with 34794 features
## It has 6 fields

sf_clc18 <- st_as_sf(clc_18)

sf_clc18_urb <- subset(sf_clc18, CODE_18 == "111" | CODE_18 == "112") %>% # CLC urban zones
  st_transform(crs = "+init=epsg:32634")

## Warning in CPL_crs_from_input(x): GDAL Message 1: +init=epsg:XXXX syntax is
## deprecated. It might return a CRS with a non-EPSG compliant axis order.

Sys.setlocale(locale = 'Serbian (Latin)')

## [1] "LC_COLLATE=Serbian (Latin)_Bosnia and Herzegovina.1250;LC_CTYPE=Serbian (Latin)_Bosnia and Herzegovina.1250;LC_MONETARY=Serbian (Latin)_Bosnia and Herzegovina.1250;LC_NUMERIC=C;LC_TIME=Serbian (Latin)_Bosnia and Herzegovina.1250"

opstine <- readOGR("D:/R_projects/Spatialization/Data/opstine/gadm36_SRB_2.shp", 
                   use_iconv=TRUE,  
                   encoding = "UTF-8")

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\R_projects\Spatialization\Data\opstine\gadm36_SRB_2.shp", layer: "gadm36_SRB_2"
## with 161 features
## It has 13 fields

sf_opstine <- st_as_sf(opstine)

toplane <- readxl::read_xls(path = "D:/R_projects/Spatialization/Data/toplane/Toplane_2015.xls") %>%
  mutate(GRAD = str_to_title(GRAD))

sf_opstine$SDG <- toplane$`Broj domacinstava prikljucenih na SDG`[match(sf_opstine$NAME_2, toplane$GRAD)] 

domacinstva <- readxl::read_xls(path = "D:/R_projects/Spatialization/Data/Domacinstva_2015.xls", sheet = "Stanovnistvo8")
lcl(loc = "C")

## [1] "LC_COLLATE=Serbian (Cyrillic)_Serbia.1251;LC_CTYPE=Serbian (Cyrillic)_Serbia.1251;LC_MONETARY=Serbian (Cyrillic)_Serbia.1251;LC_NUMERIC=C;LC_TIME=Serbian (Cyrillic)_Serbia.1251"

domacinstva <- cyr_lat(domacinstva)

## Warning in cyr_lat(domacinstva): Locale settings changed!

## [1] "LC_COLLATE=Serbian (Latin)_Serbia.1250;LC_CTYPE=Serbian (Latin)_Serbia.1250;LC_MONETARY=Serbian (Latin)_Serbia.1250;LC_NUMERIC=C;LC_TIME=Serbian (Latin)_Serbia.1250"

names(domacinstva) <- cyr_lat(names(domacinstva)) 

## Warning in cyr_lat(names(domacinstva)): Locale settings changed!

## [1] "LC_COLLATE=Serbian (Latin)_Serbia.1250;LC_CTYPE=Serbian (Latin)_Serbia.1250;LC_MONETARY=Serbian (Latin)_Serbia.1250;LC_NUMERIC=C;LC_TIME=Serbian (Latin)_Serbia.1250"

domacinstva <- domacinstva %>%
  mutate_all(~replace_na(., 0))
domacinstva$Opština[domacinstva$Opština == "Indjija"] <- "Inđija"
domacinstva$Opština[domacinstva$Opština == "LJubovija"] <- "Ljubovija"
domacinstva$Opština[domacinstva$Opština == "Mali Idjoš"] <- "Mali Iđoš"
domacinstva$Opština[domacinstva$Opština == "Savski venac"] <- "Savski Venac"
domacinstva$Opština[domacinstva$Opština == "Stari grad"] <- "Stari Grad"
domacinstva$Opština[domacinstva$Opština == "Petrovac na Mlavi"] <- "Petrovac"
domacinstva$Opština[domacinstva$Opština == "Arandjelovac"] <- "Aranđelovac"
domacinstva$Opština[domacinstva$Opština == "LJig"] <- "Ljig"
domacinstva$Opština[domacinstva$Opština == "Žitoradja"] <- "Žitorađa"
domacinstva$Opština[domacinstva$Opština == "Medvedja"] <- "Medveđa"

sf_opstine$Br_domacinstva <- domacinstva$Ukupno_domaćinstva[match(sf_opstine$NAME_2, domacinstva$Opština)]

sf_opstine %<>% 
  st_transform(crs = "+init=epsg:32634") %>% 
  mutate_all(~replace(., is.na(.), 0)) %>%
  mutate(Br_domacinstva_SDG = Br_domacinstva - SDG)

sf_clc18_urb <- st_join(sf_clc18_urb, sf_opstine, largest = TRUE) 

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

sf_clc18_urb %<>% dplyr::select(.,Br_domacinstva_SDG, NAME_2)
sf_clc18_urb[,vars] <- NA
sf_clc18_urb %<>% sf::st_transform(4326)
sf_clc18_urb.int <- st_intersection(sf_clc18_urb, sf.grid.5km) %>% 
  filter(!is.na(Br_domacinstva_SDG))

## although coordinates are longitude/latitude, st_intersection assumes that they are planar

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

source.1A4bi$sources$polygon <- sf_clc18_urb.int

sf.1A4bi <- corsum2sf_polygon(source.1A4bi, distribute = FALSE) 

## Note: Using an external vector in selections is ambiguous.
## i Use `all_of(vars)` instead of `vars` to silence this message.
## i See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.

sum_s <- sum(sf.1A4bi$Br_domacinstva_SDG)
diff.1A4bi <- data.frame(total.1A4bi - sum.1A4bi)
sf.1A4bi <- sf.1A4bi %>%
  mutate(NOx = ((diff.1A4bi$NOx/sum_s)*Br_domacinstva_SDG),
         SO2 = ((diff.1A4bi$SO2/sum_s)*Br_domacinstva_SDG),
         PM10 = ((diff.1A4bi$PM10/sum_s)*Br_domacinstva_SDG),
         PM2.5 = ((diff.1A4bi$PM2.5/sum_s)*Br_domacinstva_SDG),
         NMVOC = ((diff.1A4bi$NMVOC/sum_s)*Br_domacinstva_SDG),
         NH3 = ((diff.1A4bi$NH3/sum_s)*Br_domacinstva_SDG))
sf.1A4bi %<>% dplyr::select(vars)

p.1A4bi <- sf.grid.5km %>%
  st_join(sf.1A4bi, join = st_contains) %>% 
  group_by(ID) %>%
  summarize(NOx = sum(NOx, na.rm = TRUE),
            SO2 = sum(SO2, na.rm = TRUE),
            PM10 = sum(PM10, na.rm = TRUE),
            PM2.5 = sum(PM2.5, na.rm = TRUE),
            NMVOC = sum(NMVOC, na.rm = TRUE),
            NH3 = sum(NH3, na.rm = TRUE)) %>% 
  mutate(ID = as.numeric(ID))

## although coordinates are longitude/latitude, st_contains assumes that they are planar

## `summarise()` ungrouping output (override with `.groups` argument)

sum.p.1A4bi <- p.1A4bi %>% 
  st_drop_geometry() %>%
  dplyr::select(., vars) %>% 
  apply(., 2, sum) %>% 
  t(.) %>% 
  as.data.frame() %>%
  dplyr::mutate_if(is.numeric, round, 2)
data.frame(sum = c("spatialized", "total", "diff"), rbind(sum.p.1A4bi, total.1A4bi, data.frame(sum.p.1A4bi == total.1A4bi)-1)) %>%
  datatable(., caption = 'Table 6: Summary differences after spatialisation',
            options = list(pageLength = 5)
  )